Molecular Dynamics

Syma Khalid

Velocity Verlet algorithm

Considers both the atomic positions and velocities at the same time.

Divide the procedure into three steps:

$$ \begin{align} v(t + \delta t / 2) & = v(t) = \frac{a(t) \, \delta t}{2}, \\ r(t + \delta t) & = r(t) + v(t + \delta t / 2) \, \delta t \\ v(t + \delta t) & = v(t + \delta t / 2) + \frac{a(t + \delta t) \, \delta t}{2}. \end{align} $$

Only set of positions, forces and velocities need to be carried at any one time.

Advantages

• Simple to code
• Relatively modest memory requirements
• Velocities explictly calculated, thus temperature is reliably calculated
• Energy well conserved.

The algorithm of choice for most popular MD codes.

Temperature

Temperature specifies thermodynamic state - important concept in MD:

$$ \begin{equation} T(t) = \sum_{i=1}^N \frac{m_i v_i^2(t)}{K_B (3 N -3)} \end{equation} $$

$N$ is the number of particles and $K_B$ the Boltzmann constant.

Timesteps

How to choose $\delta t$? Must strike balance between conserving energy and efficient simulation.

High frequency motions can be problematic; get around this by "freezing out" - constraining bonds to their equilibrium values.

For biomolecules $\delta t = 2$fs is often used.

System setup

Need coordinates for the molecules to start.

Usual to get structure from x-ray crystallography or NMR (for biological systems).

Embed structure in an environment that mimics experimental conditions.

Set system up carefully, taking care to remove any "bad contacts" - eg, use energy miinimisation.

Periodic Boundaries

Only small numbers of molecules in a cell (max 150k as of 2014). Molecules on the surface experience different forces to bulk molecules, so to avoid this use periodic boundaries.

If a particle leaves the central box one of its images will enter the box through the opposite face.

The number of particles in the central box (which is actually simulated) remains constant.

System Equilibriation

The equivalent of a "well-stirred" solution in synthetic chemistry.

When doing a simulation in the NPT (Number of particles, Pressure, Temperature: these are forced to remain constant) ensemble, then we monitor the system volume (ie, the volume of the box).

Similarly, when using an NVT (Volume...) ensemble, we monitor the system pressure.

Often start by equilibriating in the NVT ensemble, then switching to the NPT ensemble.

Can study the systematic drift in the structure by monitoring the root-mean-square deviation

Timescale limitations

One of the major limitations: timescale that can be simulated.

State-of-the-art: simulations max 1-10 $\mu$s of "real time".

Problem as conformation changes can take ms, s, to hours or longer.

Additional problem: encounters between two molecules can lead to changes on diffusion timescale of ms or longer.

Enhanced Sampling Methods

Many different methods to tackle the timescale limitations. Focus on Metadynamics.

Assume system can be described by a few selected degrees of freedom, called collective variables (CV).

Sampling facilitated by addition of a history-dependent bias potential that acts on the CVs.

Calculate the location of the system in the space determined by the CVs, add a positive Gaussian to discourage the system from revisiting that location (as it has already been sampled).

How to choose CVs?

They are a function of the microscopic coordinates; not unique.

• Should distinguish between initial and final state and also describe key intermediates (but none of these are necessarilly known!).
• Should include all slow modes of the systesm.
• Should be limited in number.

Once full energy landscape sampled, the free-energy is a constant as a function of the CVs. Can construct the true free-energy by substracting the Gaussians added (big additional benefit).

How many Gaussians to add, and the width, are tuneable parameters (plenty of trial and error).

Energy minimisation

Method employed to bring system to its (local) energy minimum. To find the minimum points

• map the complete hypersurface (potential energy wrt all possible parameters) and thereby find all possible values
• want approach that keeps computational cost realistic
• tyipcally a $3N-5$ or $3N-6$ dimensional problem.

Various algorithms that work.

• No single algorithm always works
• Algorithms "go downhill", locating the nearest minimum (not necessarily the global minimum).

Two general categories:

• Those that use derivative information
• Those that don't.

Focus on one example of a derivative method.

First Derivative Method - Steepest Descent.

Function must be continuous and differentiable. First order methods use the gradients. Clearly related to the force $\vec{F} = -\nabla V$.

Imagine standing on top of a hill and looking for the steepest slope that takes you down.

• At each interation the search direction is taken as the negative gradient of the function.
• Negative goes downhill.
• How far to go? Use line search.

• Locate minimum along the specified direction.
• Bracket the minimum. Then use Bisection or similar to find the minimum.
• From this point, recalculate the search direction; must be perpendicular to the first.

Advantages:

• Quickly removes worst steric clashes
• Rapidly brings bond lengths and angles to ideal values

Disadvantages

• Inefficient in valleys and near the minimum